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SUMMARY 


This paper describes the Langley Research Center Orbital Lifetime program 
which analyzes motion of Earth-orbiting spacecraft based on perturbations due 
to atmospheric drag, solar radiation pressure, and gravitational effects of 
the Sun, the Moon, and Earth oblateness. The program can be used to generate 
time histories of the orbital elements, total lifetime and decay rates, and 
plot data. Operation Instructions, descriptions of the input and output, and 
several sample cases are included in this paper. 
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SYMBOLS 

Definition 

O 

Area of spacecraft projected perpendicular to velocity, cr.r 

O 

Acceleration of spacecraft due to atmospheric drag, cm/sec 

Geomagnetic index, dimensionless 

area of spacecraft projected perpendicular to Sun's 
direction, cm 2 

Acceleration of spacecraft due to solar radiation pressure, 
cm/sec 2 

Inertial azimuth, measured positive clockwise from north, 
deg. 

Semimajor axis of orbit, R e 

0L program FORTRAN array containing floating point input 
data 

Speed of light, cm/sec 

Spacecraft drag coefficient, dimensionless 

Eccentricity of orbit, dimensionless 

10.7 cm solar flux, 10" 22 watts/m 2 /Hz 

13-month averaged 10.7 cm solar flux, 10“ 22 watts/m 2 /Hz 

Height of apogee above Earth, km 

Height of perigee above Earth, km 

Solar flux near the Earth, dynes/cm/sec 

Inclination of orbit, deg. 
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Earth oblateness coefficients corresponding to the second, 
third, fourth, and fifth harmonics, respectively, 
dimensionless 

Spacecraft reflection coefficient, dimensionless 

OL program FORTRAN array containing integer input data 

Spacecraft mass, g 

Semilatus rectum of orbit, R e 

OL program FORTRAN array (see ref. 2) 

Earth radius, 6.378145 x 10 R cm 
Initial geocentric radius, km 

Spacecraft velocity in direction of flight, cm/sec 
Spacecraft inertial velocity, km/sec 
Inertial flight path angle, deg. 

Change in semimajor axis, R e 

Change in eccentricity, dimensionless 

Change in semilatus rectum, R e 

Change in spacecraft velocity, R p /sec 

Initial longitude, positive east of Greenwich, deg. 

Atmospheric density, g/cm J 

Standard deviation 

Initial latitude, deg. 

Right ascension of ascending node, deg. 

Argument of perigee, deg. 


INTRODUCTION 

The Langley Research Center (LaRC) Orbital Lifetime (OL) program analyzes 
the long-term motion of Earth-orbiting spacecraft for altitudes up to 
2500 km. It models perturbations to the orbit caused by atmospheric drag, 
solar radiation pressure, and gravitational effects due to the Sun, the Moon, 
and the oblateness of the Earth. Available output includes time histories of 
the orbital elements, total lifetime and decay rates, and plot data. 
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The OL program was originally developed at the Massachusetts Institute of 
Technology Lincoln Laboratory (ref. 1) and was modified by the Rand 
Corporation (ref. 2). NASA LaRC obtained the program from Rand and further 
changes have been made by LaRC personnel. These changes Include correction of 
program errors, updating to later versions of FORTRAN, and replacing the 
atmospheric density model. This user's guide and its references constitute 
all of the documentation available for the OL program. The theory in 
reference 2 is still reliable, but most information in that document specific 
to the program code itself is no longer accurate because the code has been 
changed. Although this paper accounts for the changes that were made after 
NASA LaRC obtained the program, the program still contains some calculations 
which were originally in the code but were never documented. 

The purpose of this paper is to provide a user's guide that is consistent 
with the current version of the OL program. This paper gives, without going 
into very much mathematical detail, an overview of the method of computation 
and the modeling approach and underlying assumptions for each perturbation. 

As mentioned above, the theory behind the program (i.e. the equations of 
motion of the orbital elements and derivations of expressions for the 
perturbing accelerations) appears in reference 2. The atmospheric density 
model is documented in references 3-5. This paper includes instructions for 
setting up the input and running the program, as well as sample cases which 
illustrate many of the options available. The OL program is available in 
FORTRAN 7/ on VAX and Prime computers and in FORTRAN V on CDC Cyber 
machines. The instructions and sample cases in this paper refer to VAX 
computers and the VAX/VMS operating system. 

This paper is intended to provide all of the documentation needed for the 
general user to run the OL program with confidence. 


PROGRAM DESCRIPTION 

The OL program calculates changes in spacecraft orbits, given such input 
quantities as spacecraft physical characteristics, initial orbit parameters, 
and launch date. The calculations begin with the assumption that the 
spacecraft is in the initial orbit specified, so launch date is actually 
treated as initial time in orbit. The orbit is allowed to decay due to 
environmental perturbations, and no reboost operations can be taken into 
account. (There is a seldom-used option that allows for a aV kick at perigee; 
th's is not reboost.) For purposes of this program, only long-term changes 
are considered; short-period variations (e.g. within a single orbit) are 
assumed to average out and are not included in the calculations. Therefore, 
the OL program cannot be used to investigate short-term behavior of spacecraft 
flight. The perturbations to the orbit are calculated to first order only 
(except for Earth oblateness, which can include a second-order term in the 
gravity model; see Earth Oblateness subsection below). 

The spacecraft orbit is represented in the program by the orbital 
elements p (semilatus rectum), e (eccentricity), i (inclination), « (argument 
of perigee), and ft (right ascension of the ascending node). Epoch (time of 
perigee passage) is not included because the program is not concerned with the 
exact position of the spacecraft in its orbit at any given time. Changes in 


the orbital elements are calculated based on solutions to their equations of 
motion expressed in terms of the perturbing accelerations and the orbital 
elements themselves. The equations, which are of the form 

d (element) * f{p, e, i , u>, Q, perturbing accelerations), 
dt 

are integrated over a single orbit, holding constant the orbital elements 
which appear on the right-hand side of the equation. Each perturbation is 
handled separately (see descriptions in the following section). The 
integrations are analytic for all of the perturbations except drag, for which 
an approximation technique is used. The integration results are used to 
update the orbital elements at user-specified intervals (input as number of 
orbits per calculation), assuming that their rates of change remain constant 
over the computation interval. The updated orbital elements then describe a 
new spacecraft orbit for which the equations c* motion are again integrated. 
The cycle continues until a user-specified maximum orbit number is reached or 
the spacecraft reaches the Earth. The condition for reaching the Earth is 
that the height of perigee be below 64 km; the atmosphere is so dense at that 
altitude that the spacecraft would fall from there to the ground in a matter 
of minutes. 

Appendix A contains a brief description of each routine in the program. 


PERTURBATIONS 

A brief explanation of each perturbation is provided here. For more 
mathematical detail, see reference 2. 

Solar Radiation Pressure 


In the OL program, solar radiation pressure is represented by an 
acceleration of the spacecraft in a radial direction away from the Sun: 

A rp * k (I/c) (A r /M) 


where 

A RP = acceleration of spacecraft due to solar radiation pressure 

k = spacecraft reflection coefficient; 0 _< k _< 4/3 
k * 0: transparent 

k * 1: perfectly absorbing 

k = 4/3: flat, specularly reflecting 

I * solar flux near the Earth 

c 3 speed of light 

A r = area projected perpendicular to Sun's direction 
M * spacecraft mass 

The solar flux, I, is assumed to vary only with the Earth's distance from 
the Sun (the Earth-Sun distance varies because the Earth's orbit is 
elliptical). The program uses I/c - 4.5 x 10 s dynes/cnr at an Earth-Sun 
distance of one astronomical unit. 
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The effect of the Earth ^ shadow Is taken into account in calculating the 
radiation pressure by excluding from the orbital integration that portion of 
the orbit, if any, which lies in the Earth's shadow. The true anomalies on 
entering and leaving the shadow are printed in the output, and there is an 
option to calculate and print the amount of time the spacecraft spends in the 
shadow. 

Subroutine RADPR calculates the solar radiation pressure effects and 
subroutine SHADOW calculates the true anomalies on entering and leaving the 
Earth's shadow. 

Earth Oblateness 


The OL program includes the gravitational effects of the oblateness of 
the Earth through an expansion of the gravitational potential in spherical 
harmonics. Zonal harmonics are included through Jc. It is assumed that the 
Earth's mass distribution has azimuthal symmetry (i.e.that it is independent 
of longitude), so no tesseral harmonics are used. 

As mentioned above, there is an option to include a second-order ^ 
correction in the Earth oblateness calculation. This correction is the J? 
term, and this option was included because it contributes as much to the K 
gravitational potential as the J3, J4, and J5 terms do. 

It should be noted that the changes in the semimajor axis, a, due to 
Earth oblateness should be very small compared to the changes due to drag. 
However, the Ap and Ae originally calculated in the program gave rise to long- 
term oscillatory variation in a (a = p/(l-e^)) with Aa per orbit of the same 
order of magnitude as Aa due to drag. The reason for this error has not been 
found. Therefore, Ap and Ae were set equal to zero in the program to yield 
more accurate results. 

The Earth oblateness calculations appear in subroutine EARTH. 

Moon's Gravity 

The OL program accounts for the gravitational attraction both between the 
Moon and the spacecraft and between the Moon and the Earth. That is, it 
calculates the acceleration of the spacecraft due to the Moon's gravity 
relative to the acceleration of the Earth due to the Moon's gravity. The 
calculations account for distance variations due to the Moon's orbital motion. 
Terms linear in the ratio of Earth-spacecraft distance to Earth-Moon distance 
are included, but higher powers of that ratio are not. 

Subroutine MOON contains the computation for the Moon's gravity. 

Sun's Gravity 

The gravitational effects of the Sun on the orbit are handled in the OL 
program in a manner analagous to the method used for the Moon's gravity 
calculations. 

Subroutine SUN contains the calculations for the Sun's gravitational 
effects. 
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Atmospheric Drag 

The OL program calculates the acceleration of the spacecraft due to 
atmospheric drag In a direction opposite to the spacecraft velocity according 
to: 

A D r * U/2) (C 0 A 0 /M) P V 2 

where 

Aq R * acceleration of spacecraft due to atmospheric drag 
Aq * area projected perpendicular to the velocity, averaged over 
one orbit 

C D = spacecraft drag coefficient 
M = mass of spacecraft 
p = atmospheric density 
V = velocity of spacecraft 

The equations of motion in which this expression appears cannot be 
integrated analytically because the atmospheric density, p, appears inside the 
integrand, and it varies throughout the orbit in a complicated way (see below 
for a description of the atmospheric density model used). Therefore, the 
integration is approximate; six 9-point Gaussian quadratures are used, and 
there is an option to change the number of quadratures. An oblate Earth is 
used for the drag calculation, i.e. the height at which the atmospheric 
density is calculated is the height above an oblate rather than a spherical 
Earth. In addition, there is an option to use a rotating atmosphere. 

The atmospheric density models used in the OL program's drag calculations 
are the U.S. Standard Atmosphere 1976 (ref. 3) for altitudes below 90 km and 
the Jacchia 1970 atmospheric density model (refs. 4 and 5) for altitudes above 
90 km. The Jacchia 1970 model is recommended by NASA for use in space vehicle 
design and development studies (ref. 6). These models are new to the LaRC 
program and replace the models used In the Rand version of the program 
(ref. 2). 

The Jacchia 1970 atmospheric density model of the neutral atmosphere is 
valid from 90 kr to 2500 km. It includes variations with altitude, solar 
activity (I.e. with the 11-year solar cycle), and geomagnetic activity. It 
also includes the semiannual variation (characterized by maxima in April and 
October and minima in January and July), the diurnal variation (density is 
higher on the light side of the Earth than on the dark side), and seasonal - 
latitudinal variations in helium density and total density. The OL program 
calculates the seasonal -latitudinal variations based on the Jacchia 1971 
atmospheric density model (ref. 7). 

To account for solar activity and geomagnetic activity, the Jacchia 1970 
model uses the 13-month smoothed 10.7 cm solar flux, P.q ,, and the geomagnetic 
index, A . Density increases as either of these increases, but 7 has 
a greater effect. Since the OL program Is used primarily to predict 'lifetimes 
for_future flights of spacecraft, it uses predictions of monthly values 
of 7 and A . These predictions are provided by the Atmospheric Sciences 
Division of the NASA Marshall Space Flight Center (MFSC) and are updated 
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bi-monthly. The OL program reads these predictions from a flux input file, 
and both nominal (50th percentile) and +2 <j (97.7th percentile) predictions are 
available for use. A listing of the flux input file containing the MSFC 
predictions from the June 1985 update is provided in Appendix B. The earliest 
date for which solar flux data exists in the flux file is the earliest launch 
date for which the program can be run (January 1984 for the file in Appendix 
B). The date of the last flux data point (November 1998 for the file in 
Appendix B) presents no such restriction because for dates beyond the end of 
the flux data, the 11-year solar cycle is assumed to repeat as many times as 
necessary. 

The atmospheric drag calculations appear in the program in subroutine 
DRAG1, and atmospheric densities are calculated in function RHA. 


HOW TO RUN THE OL PROGRAM 

These instructions for compiling and running the OL program apply to VAX 
computers with the VAX/VMS operating system. The sample cases below were run 
on a VAX 11/785 with the VMS 4.1 operating system. 

Compiling the OL program is straightforward because there are no 
libraries or special options required. If the FORTRAN source code is in the 
file OL.FOR, the following sequence of commands compiles, links, and runs the 
program: 

FOR OL Compiles the highest-numbered version of OL.FOR and creates 

a new version of OL.ORJ 

LINK OL Links the highest-numbered version of OL.OBJ and creates a 

new version of OL.EXF 

RUN OL Runs the highest-numbered version of OL.EXF 

There are two input files required to run the program: 

1. Case-specfic input data file: This file, created by the user, 

contains information such as initial orbit parameters and case 
flags. Its format is described in detail in the Input Format section 
below. The input datj file may have any name of the form 
fi1ename.DAT, which the user specifies at the beginning of each run 
(see example below). Two cases using the same input file may not be 
run simultaneously because this file is open during the entire run. 

2. Flux data file: This is a file named FLUX.DAT which contains 

predictions of solar flux and geomagnetic index needed for the 
atmospheric density model. Appendix B contains a copy of the flux 
data file used in the sample cases below. The flux file is read at 
the start of each run and is closed immediately after it is read. If 
the flux file is already open when the program tries to read it, the 
program tries again until it becomes available. Therefore, multiple 
cases may be run simultaneously using the same flux file. 

The two input files must exist in the directory under which the program is 
run. 
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The OL program generates two output files: one contains the output 
listing and one contains data for plotting (see sample cases below). The user 
specifies names for the output listing and plot files at the beginning of the 
run (see example below). 

When the program Is run Interactively, It prompts the user for relevant 
file names at the terminal. The following Is a sample terminal session which 
runs the program file OL.EXE, using Input file OLIN.OAT and flux file FLUX. OAT 
(required name). The run generates output listing file 0L0UT.DAT and plot data 
file OLPLT.OAT. Lower case letters Indicate user Input. $ Is the system prompt. 

$ run ol 

ENTER "INPUT FILE" NAME, 
olln 

ASSIGNING FLUX FILE. 

ENTER "OUTPUT FILE" NAME, 
olout 

ENTER "PLOT FILE" NAME, 
ol pi t 
5000 
$ 

The three file name prompts ap;*ar at the beginning of program execution. The 
number 5000 appears when the program terminates normally. If the specified 
input file does not exist, execution ends and 

FILE ERROR 
$ 

appears in place of the remaining prompts. If the flux file does not exist, 

EOF FLUXIN 

appears In place of 5000. If any other error occurs during execution, the 
appropriate system error message appears. 

It Is recommended that the OL program be run In batch mode because 
execution times are typically too long for Interactive execution of the 
program to be convenient (see sample cases below), especially for multiple 
cases. To run in batch mode, a command file must be created and submitted. 

The following command file, named OL.COM, sets up the case in the above 
interactive example to execute under the directory LHO.OL: 
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$SET DEFAULT [LHO.OL] 

SR UN OL 

OLIN 

OLOUT 

OLPLT 

SEXIT 

The command 
SUBMIT OL 

submits the highest numbered version of the file OL.COM for a batch run. 


INPUT FORMAT 


There are four types of data in the OL program input file. They are, in 
order of appearance: 

1. Initial line which specifies run number and the numbers of lines each 
of character, integer, and real data that follow; 

2. Lines containing character data to be used as an output heading; 

3. Lines containing integer data such as output flags; and 

4. Lines containing real data such as orbit specifications. 

A sample input file for one of the five sample cases described below can be 
seen in Table 1. 


Initial Line 


All values are integer and must be right justified in their fields. 


Column #'s Default 


Definition 


1-6 

25-30 

31-36 

55-60 


0 Run number. Used only as a label on the output. 

0 Number of lines containing Integer data. 

.0 Number of lines containing real data. 

0 Number of lines containing character data. 


Columns 7-24, 37-54, and 61-72 are ignored. 

Character Data 

Up to two 72-character lines of text may be input to be used as an output 
heading. There is no default heading. 
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Integer Data 

All Input integers are assigned to elements of an array named K in the 
program. The array index and the value of the input integer must be listed 
consecutively, with six columns for each, right justified. Read format is 
1216, so up to six such pairs may be included in each line. The order of 
input is arbitrary; e.g. K(10) may precede K(2). Blank pairs of fields are 
ignored. 

The most commonly used elements are: 

K Index Default Definition 


1 


2 

5 


10 


50 


54 


2 Maximum number of orbits in case. Case ends 

when K(l) orbits are reached or when spacecraft 
reaches Earth, whichever occurs first. 

K(l) = 999999 gives whole lifetime for most 
cases. 

1 Number of orbits per calculation. 

6 Number of Gaussian quadratures in drag 

calculation. Number of calculation points per 
orbit = K(5)*9. 

1 Print interval. Print every <(10) orbits. 

Actual print interval is least common multiple 
of K ( 2 ) and K( 10) . 

0 Orbit specification option. 

K(50) = -1: Input C( 1108), C(1110), C(llll), 

C(1112), C ( 1 1 1 3 ) , C ( 1 1 14 ) 

K(50) = 0: Input C(537), C(538), C(539), 

C(540), C ( 54 1 ) 

K ( 50 ) = 1: Input C(539), C(541), C(llll), 

C( 1121 ), C( 1122) 

The C array contains real input data and is 
described below. 

0 Solar flux flag. 

K(54) = 0: Nominal solar flux and geomagnetic 

index predictions used. 


K(54) = 1: +2a predictions used. 
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Less often used elements of K are: 

K Index ~ ault Definition 

14 0 Orbit number at which print interval changes 

to K(15). 

15 0 Print interval after orbit number K(14). 

71 0 Gravity model flag for Earth oblateness 

calculation. 

K ( 71 ) = 1: Jp term included. 

o 

K ( 71 ) = 0: jp tern not included. 

82 0 Shadow flag. 

K(82) = 1: Time in shadow per orbit printed. 
K(82) = 0: Tine in shadow not printed. 

111-118 0 Flags to delete perturbations. Setting any one 

of K(lll) - K ( 1 1 R ) equal to one of the numbers 
below deletes the effect of the perturbation 
associated with that number: 


Solar radiation pressure = 5 
Earth oblateness = 7 
Moon's gravity = 8 
Sun's gravity = 9 
Drag = 6 


Using this option does not save any computation 
time; changes due to the specified 
perturbations are calculated in the usual way 
and then set equal to 0, 

The following elements of K generate more detailed output and are rarely 

used: 

K Index Default Definition 

18 0 Flag for number of lines of printout per print 

interval (K(10) or K ( 1 5 ) ) . 

K(18) = 0 or 1: One line 

K ( 18 ) * 2: Two lines 

K ( 1 8 ) > 3: K ( 18) lines 

K(18) 7 -3: -K(18) lines 

See ref. 2 for more details. 


i • 
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20 


21 

22 

24 


25-29 


0 Print interval for long output, which can 

include contributions of each perturbation to 
changes of the orbital elements; see K(24) 
below. Actual long output print Interval 3 
least common multiple of K(2), K(10), and 
K(20). This long output Is In addition to 
normal output controlled by K(10). 

0 Orbit number to begin long output. 

0 Orbit number to end long output. 

0 Flag for number of lines of long output. 

K(24) * 0: seven lines 

K(24) > 0: K(24) lines 

The first two lines contain parameters 
describing the orbit. The next five lines list 
the contributions of each perturbation to 
changes in the orbital elements. There is no 
useful Information beyond the seventh line. 

0 Row of Q array to be included in long 

output. See reference 2. 


Real Data 


All Input real numbers are assigned to elements of an array named C in 
the program. As with integers, the array index and value of the element are 
listed in pairs. There are 4 columns for the array index, which is intege' 
and right justified, and 14 columns for the value of the element, which is 
real. Blanks are read as zeroes. Read format is 4(14, E14.7), so up to four 
pairs will fit on a line, and the order of input is arbitrary, e.g. C(llll) 
may p ecede C(539). Blank pairs of fields are ignored. 

The most commonly used elements are: 

Spacecraft Characteristics 

C Index Units Default Definition 

531 cm 2 /g 0.0 kA R /M, where 

k = spacecraft reflection coefficient 
k = 0: transparent 

k * 1: perfectly absorbing 

k * 4/3: flat, specularly 

reflecting 

A R * area projected perpendicular 
to direction of Sun, cnr 
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532 cm 2 /g 0.0 


M = total mass of spacecraft, a 

If C ( 531 ) < 10" fi , radiation pressure 
calculation is omitted. 


C 0 A n /2M, where 


C D * spacecraft drag coefficient 
Ap = area projected perpendicular to 
velocity, cm 2 

M = total mass of spacecraft, g 


If C( 532 ) < 10“ fi , drag calculation 
is omitted. 


Orbit Specification 

There are three ways to specify the initial orbit, given by K(50) = 0, 1, 
or -1. There is some redundancy among the three options (* indicates element 
common to two options). The first two options, K(50) = 0 or 1, assume that 
the initial condition is at perigee. For the third option, K(50) = -1, the 
p.-ogram calculates the time since perigee and backs up the initial time to 
correspond to perigee. If either of the two options K(50) = ±1 is selected, 
the program calculates the orbital elements corresponding to K(50) = 0; see 
Appendix C. Note that the maximum allowed altitude is 2500 km. 


For 

K(50) = f 

): 


C Index 

Units 

Default 

Definition 

537 

R e 

0.0 

p = semilatus rectum of orbit 

538 

— 

0.0 

e » eccentricity of orbit 

539* 

deg 

0.0 

«i) = argument of perigee 

540 

deg 

0.0 

n = right ascension of ascending node 

541* 

deg 

0.0 

i * inclination of orbit 


For K(50) 

= 1: 


C Index 

Units 

Default 

Definition 

539* 

deg 

0.0 

:) = argument of perigee 

541* 

deg 

0.0 

i = inclination of orbit 

1111* 

deg 

0.0 

X = longitude 

1121 

km 

0.0 

hp = height of perigee above Earth 

1122 

km 

0.0 

h a = height of apogee above Earth 
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For K(50) « -1: 

C Index Units Default 


Definition 


1108 

deg 

0.0 

Az 

* inertial azimuth 

1110 

deg 

0.0 

$ 

= latitude 

mi* 

deg 

0.0 

X 

* longitude 

1112 

km 

0.0 

R i 

= geocentric radius 

1113 

km/sec 

o.n 

v i 

= initial velocity 

1114 

deg 

0.0 

y 

= flight path angle 


Launch Date 


Launch date (initial time) can be input either as Julian date (C ( 542 ) ) or 
as calendar date (C(1115)-C(1120) ): 


C Index 

Units 

Default 


Definition 

542 

days 

0.0 


Julian date 

1115 

years AD 

0.0 


Year 

1116 

months 

0.0 


Month 

1117 

days 

0.0 


Day 

1118 

hours 

0.0 


Hour 

1119 

minutes 

0.0 


Minute 

1120 

sec 

0.0 


Second 

Seldom used elements 

of C are: 



C Index 

Units 

Default 


Definition 

564 

-- 

0.0 

Rotating atmosphere flag. 




C(564) = 
C(564) - 

1.0: Uses rotating atmosphere 

0.0: Uses ron-retating 

atmosphere 

565 

R e 

.0004 

Change in semimajor axis of orbit (t&) < 
which to reduce computation and print 


intervals to one. When a a >C ( 565 ) from 
one orbit to the next, K(2) and K(10) are 
set equal to one. 
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566 — 0.0 Flag for computing gamma in drag 

calculation (see ref. 2. for definition 
of gamma). 

C ( 566 ) = 1.0: Gamma computed 

C(566) = 0.0: Gamma = 1.0 

When gamma is computed, it should appear 
in the integrand in the orbital 
integrations. However, in the program, it 
appears outside the integral as a 
multiplier, so choosing to have gamma 
computed does not necessarily give better 
results than letting gamma = 1. 

570 R e /sec 0.0 aV kick at perigee. 


SAMPLE CASES 

Tables 1, 2, and 4-11 contain input and output files for each of five 

sample cases. Table 3 is a listing of the plot data file for case 1, and fig. 

1 is the corresponding lifetime plot. The first four cases represent the NASA 
Space Station Reference Configuration (ref. 8) in a 500 km (270 nmi ) circular 
orbit. Case 5 represents a spherical satellite in a 600 km by 900 km polar 

orbit. Each case illustrates several input and output options. 

Case 1 


Table 1 shows the input file for sample case 1. Integer data for the K 
array appear in the fourth line: 


K( 1 ) = 999999 

= Maximum orbit number: gives entire lifetime. 

o 

*-H 

It 

CM 

Calculates every 10 orbits. 

K( 10) = 100: 

Prints every 100 orbits. 

K ( 7 1 ) = 1: 

p 

Includes Jp term in Earth oblateness calculation. 

K( 50) = 1: 

Selects orbital input option where C(539), 

C ( 541 ) , C(llll), C( 1121 ) , and C(1122) must be input 

K(54) = 1: 

+ 2a flux predictions used. 


Real data for the C array appear in the fifth through seventh lines: 
Spacecraft: 

C ( 531 ) = 0.1543 = kA^/M for reference configuration 
C ( 532 ) = 0.1365 = CqA^/ 2M for reference configuration 
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Orbit: 

C (539) = 0.0 = argument of perigee 
C ( 541 ) = 28.5 = Inclination 
C(llll) = 0.0 = longitude 

C ( 1121 ) = 500.0393 = height of perigee above Earth 

C ( 1122 ) = 500.0393 = height of apogee above Earth 

This is a circular, 500 km orbit, inclined 28.5° to the equator. 

Launch date: 

C (1116) - 6.0 = month 
C ( 1 1 1 7 ) = 1.0 = day 
C ( 1 120 ) = 1.0 = second 
C ( 1 1 15 ) = 1991.0 = year 

Launch is 1 second after midnight June 1, 1991. This launch date 

corresponds to flight during a maximum in solar activity and thus, a 

maximum in atmospheric density (see Appendix B). 

Table 2 is a listing of the output file generated for case 1. The first 

four lines on the first page repeat the K and C values and array indices that 

were input. The next printed line lists the Julian and calendar dates (GMT is 

Greenwich mean time). Next is the heading from the input (this is the only 

place where the heading appears). Following the heading is a line which 
repeats information from the first line of the input file: 

NUR = run number 

NK = number of lines containing K array data 

NC = number of lines containing C array data 

NH = number of lines containing character data 

The final three lines on the first page of output display initial values of 

some variables and some constants: 

RAM = kA R /M, cm 2 /g 

DAM = C d A d / 2M, cm 2 /g 

XP = initial semilatus rectum of orbit, R e 
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XE » Initial eccentricity of orbit (Note that although height of 

perigee and apogee are equal In the 1- out, XE is not identically 
zero. XE is set equal to 10"° in the program to avoid 
problems associated with divisio' by ;ero.) 

XOMEGD * initial argument of perigee, dev,. 

XASCND * initial right ascension of ascending node, deg. 


XOINCD = initial inclination of orbit, deg. 

XDAY = initial number of days since midnigh J , 12/30/57, days 

PSID = angle between equatorial and e-lipti. planes, deg. 

DH = altitude below which spacecraft is considered to have reached the 

Earth, 10^ m 

O O 

GM = universal gravitational constant times Earth mass, cm /sec 
RE = Earth radius, cm 

OS = angular velocity of Earth about Sun, rad/day 

O 

RAOP = solar radiation pressure at Earth s surface, dyne/cm 


The second page of printout begins a tab! t which lists various orbit 
parameters as they evolve: 


J = orbit number 

P * semilatus rectum of orbit, R g 

E = eccentricity of orbit 

OMEGAD = argument of perigee, deg. 

ASCND * right ascension of ascending node, deg. 

OINCD » inclination of orbit, deg 

M.J.D. = modified Julian date, days 

A » semimajor axis of orbit, R e 

HTPER = height of perigee above Earth, R e 

DAY = number of days since launch, day r 

ALT. = height of perigee above Earth, nmi 

V1D = true anomaly on entering Earth's shadow, deg. 

V2D = true anomaly on leaving Earth's shadow ueg. (If V1D = V2D 

= -777.7777, the spacecraft does not .inter the Earth's shadow.) 


After printing values for the in.tial orbit (whose orbit number is 0), 
the program prints information about solar flux. It then continues 
calculating and printing at the specified intervals until the next calculation 
would either bring the spacecra f c below 63 km or change the semimajor axis by 
more than 0.0004 R (default value for C ( 565 ) ) . If either of these conditions 
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Is met, the calculation and print intervals are reduced to one and the 
calculations resume. For case 1, this occurs after orbit 3560; see Table 2. 
When the spacecraft reaches the Earth, the program prints out the final orbit 
number, "SATELLITE IS DOWN," final modified Julian date, semimajor axis, 
height of perigee, and end of run statement. It also prints the total 
lifetime in days (230.5 for this case) and linear decay rate in nautical miles 
per day from initial perigee altitude to 200 nautical miles. 

Table 3 is a listing of the plot data file generated for case 1. The 
first number is the number of points the file contains. Each following line 
consists of time since launch in days and perigee altitude in nautical miles; 
these are the same numbers that appear in the output listing under "DAY" and 
"ALT." The data in Table 3 are shown olotted in figure 1. The plot itself 
was generated using a separate plotting program. 

The run time for case 1 was 14 minutes and 49 seconds (CPU time on the 
VAX 11/785). 

Case 2 


Table 4 shows the input file for sample Case 2. This is the same 
spacecraft, orbit, and launch date as case 1, with some different options. 

The orbit is specified according to K(50) = 0 rather than 1, so 

C( 537 ) - C ( 541 ) are given. Launch date is input as Julian date, C(542), 

rather than calendar date. Other differences are: 

K(82) = 1: Prints time in shadow per orbit. 

K(5) = 5: Uses five Gaussian quadratures in drag calculation. 

K(14) = 2000 and K ( 15 ) * 50: After orbit 2000 prints every 50 orbits. 

Table 5 shows the output generated for case 2. The time in shadow is now 
part of the output, and the print interval is reduced to 50 at orbit 2000. 

The total lifetime is 232.3 days, essentially the same as for case 1; the 
Slight difference is due to using five Gaussian quadratures (K(5)) for drag, 
rather than the six used by default for case 1. Another effect of using five 
Gaussian quadratures is a reduction in run time: CPU time for case 2 was 

12 minutes and 26 seconds. 

Case 3 


Table 6 shows the input file for sample case 3, which has the same 
spacecraft and orbit as cases 1 and 2, but a different launch date, January 1, 
1997. As can be seen in Appendix B, the late 1990 's are characterized by 
lower solar activity than the early 1 990 ' s , when the launch date for cases 1 
and 2 occurs. Because of this, atmospheric density is lower and the total 
orbital lifetime will be longer. Therefore, to reduce run time and amount of 
printout, the calculation interval K(2) was increased to 40 and the print 
interval K( 1 0) was increased to 200. Also, the number of Gaussian quadratures 
(K(5)) is five as in case 2. 

Table 7 shows the output listing for case 3. The total lifetime is 1408 
days, longer than cases 1 and 2, as expected. 

Run time was 19 minutes and 5 seconds. 
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Case 4 


Table 8 Is the input file for case 4. It represents the same spacecraft, 
orbit, and launch date as case 3, but uses nominal solar flux predictions 
(K( 54) = 0) rather than the +2a flux used in case 3. Nominal flux predictions 
represent an average case; +2a values represent worst case. Therefore, the 
atmospheric density is lower than in case 3, and the lifetime is longer (1950 
days), as can be seen in Table 9. Case 4 uses an additional output option: 
K(20) = 200 and K(22) = 1000, which generates long output every 200 orbits 
through orbit 1000. The long output here consists of seven lines, the first 
of which is identical to the normal output and corresponds to the printed 
column headings. The second line is not labeled; there are twelve numbers 
listed, and they are (left to right): 

O 

Distance to perigee from center of Earth, 10 cm 

Change in log of period for one orbit 

Period, days 

Zero 

Zero 

True longitude of sun, deg. 

Six Zeroes 


Each of the last five long output lines lists the changes over one orbit 
in eight parameters of the orbit due to a specific perturbation: 


RADPR = solar radiation pressure 
DRAG = atmospheric drag 
EARTH = Earth obloteness 
MOON = Moon's gravity 
SUN = Sun's gravity 


The eight changes are: 


D(P) 

D(E) 

D(OMEGAD) 

D(ASCND) 

D(OINCD) 

D(A) 

D ( LOG ( PER I OD ) ) 
D(HTPER) 


change in semilatus rectum, R 0 

change in eccentricity 

change in argument of perigee, deg. 

change in right ascension of ascending node, deg. 

change in inclination, deg. 

change in semimajor axis, R e 

change in log of period 

change in height of perigee above Earth, R 0 


Run time was 25 minutes and 56 seconds. 


Case 5 


The input file for sample case 5 is shown in Table 10. The spacecraft 
and orbit are different from the first four cases: CqAq/ 2M = kA^/M here, 

which might represent a spherical spacecraft (equal projected areas for drag 
and solar radiation pressure) with Cg = 2.0 and k = 1.0. The orbit is 
elMptical, 600 km by 900 km, and polar; i.e. the inclination is 90°. 2 Launch 
date is February 13, 1985. +2o solar flux is used (K(54) = 1), and Jp terms 

are not included in the Earth oblateness calculations ( K ( 7 1 ) = 0). Tne 
maximum orbit number ( K ( 1 ) ) is 20,000. 
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Table 11 shows the output file for case 5. From the DAY and ALT columns. 
It Is clear that this Is a very slowly decaying orbit— the perigee drops by 
less than a nautical mile In almost 4 years. This Is because the atmospheric 
density Is extremely low at this orbit's altitudes. 

Note that for some orbits, V1D and V2D are -777.7777. This Indicates 
that the spacecraft never enters the Earth's shadow during these orbits. 

As specified by the input, the run terminates after 20,000 orbits, before 
the spacecraft reaches the Earth. Therefore, lifetime and decay rate are not 
printed. A plot file is still generated. 

Run time for case 5 was 22 minutes and 32 seconds. 


COTWENTS ON PROGRAM OPERATION 


In the discussion of the sample cases, the consequences of making certain 
input choices were pointed out. A few are repeated here along with other 
information to keep in mind when running the OL program. 


The OL program takes a relatively long time to run, and thus should be 
run in batch mode rather than interactively. The run times for the sample 
cases were on the order of about 20 CPU minutes. Actual turnaround times 
averaged about 1 hour when the computer had few users. The main contributor 
to the long execution time is the integration scheme in the atmospheric 
density model (function RHA) used for the drag calculations (subroutine 
DRAG1). Therefore, one way to reduce run time is to reduce K(5), the number 
of Gaussian quadratures used in the drag calculation, as was done in sample 
case 2. Case 1 used the default six quadratures; case 2 used five quadratures 
for the same orbit and the run time was reduced by about a sixth. Some 
accuracy was lost, but the difference in total lifetimes for the two cases was 
less than 1 percent. It should be noted that the orbit was circular for these 
cases and an elliptical orbit would be more sensitive to a lower K(5). 

Another way to control run time is through the choice of the calculation 
interval K(2), since run time is approximately proportional to K(2). Its 
choice should be governed by expect J lifetime, which depends primarily on 
these factors: 

C n A n /2M: Drag increases with C n A n /2M, so the higher its value, the 

shorter the lifetime. 


Altitude: Atmospheric density decreases roughly exponentially 

as altitude Increases, so lifetime increases with 
increasing initial altitude. For very high altitudes, the 
0L program calculations may not be sensitive enough for the 
orbit to decay at all. 

Launch date: Atmospheric density varies by as much as an order of 

magnitude throughout the 11-year solar cycle. As solar 
flux Increases, density Increases, so a launch during high 
solar activity results In a shorter lifetime than a launch 
during low solar activity. By current predictions, solar 
activity will be high In the early 1990's and low In the 
late 1990's (see Appendix B for solar flux predictions 
through November 1998). 
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Nominal vs. +2o solar flux: Since nominal predicted flux values 

represent the 50th percentile and +Za values represent the 
97.7th percentile, +2o values give shorter lifetimes than 
nominal values (see sample cases 3 and 4). 

The following factors do not strongly affect total lifetime: kA R /M, 

exact day and time of launch, inclusion of J? term in Earth's gravity model, 
and initial values of X, u>, and n. 

As mentioned in the input format section, the effects of individual 
perturbations may be deleted by using K( 111 ) - K(118). Using these will not 
eliminate any calculations, so they cannot be used to reduce run time. 

However, the solar radiation pressure and drag calculations can be eliminated 
by setting kA R /M (C ( 531 ) ) and C|)A 0 /2M (C(532)), respectively, to values less 
than 10" B g/cm. 

For the drag and radiation pressure calculations, the same CqAq/ 2M and 
kA r /M values are used throughout the entire lifetime and, in particular, 
throughout each orbit. Therefore, if the projected area of the spacecraft is 
expected to change, the area input to the program must be an average area. 
Sample cases 1-4 represent the NASA Space Station Reference Configuration 
(ref. 8), which has solar arrays which track the Sun. The solar array 
projected areas used were averaged over an orbit. 

The OL program will not run correctly for dates earlier than the first 
date in the flux data file (for the flux data listed in Appendix B, that date 
is January 15, 1984), To run the program for earlier dates, the flux file 
must be changed to include those dates. Dates later than the end of the flux 
data pose no problem because the solar cycle is assumed to repeat as 
necessary. 

Multiple cases may be executed in a single run by combining input for the 
cases consecutively in one input file. This option is not recommended because 
the run time is the sum of the individual cases run times; it is faster to 
submit th cases separately and let them run concurrently. 


CONCLUDING REMARKS 

This paper described the LaRC Orbital Lifetime program, which is used to 
predict orbital lifetimes and decay rates of Earth-orbiting spacecraft. The 
input and output were described in detail, and several sample cases were 
included. Further information can be obtained from the Spacecraft Analysis 
Branch, Space Systems Division, NASA LaRC. 
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TABLE II. - SAMPLE CASE 1 OUTPUT FILE 
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APPENDIX A. OL PROGRAM ROUTINES 

The FORTRAN program units In the OL program are: 

PROGRAM LIFTIME: Main program. Opens files and calls the subroutines which 

contain the perturbation calculations. 

SUBROUTINE SETUP: Initializes constants and sets defaults. 

FUNCTION TLONG: Calculates true longitude of the Sun (used in Earth 

shadow calculations). 

BLOCK DATA: Contains altitudes and densities from 1976 Standard 

Atmosphere below 90 km. 

SUBROUTINE INPUT: Reads user's input file containing case-specific data and 

sets up initial values of case-specific program variables. Writes initial 
information to output file. 

SUBROUTINE OUTPUT: Writes output file and plot file. 

SUBROUTINE SHADOW: Calculates true anomalies on entering and leaving Earth’s 

shadow (used in solar radiation pressure caTculations). 

SUBROUTINE RAOPR: Calculates changes in orbital elements due to solar 

radiation pressure. 

FUNCTION TDIST: Calculates distance from Earth to Sun (used in solar 

radiation pressure calculations). 

SUBROUTINE EARTH: Calculates changes in orbital elements due to gravitational 

effects of Earth oblateness. 

SUBROUTINE MOON: Calculates changes in orbital elements due to gravitational 

effects of the Moon. 

SUBROUTINE SUN: Calculates changes in orbital elements due to gravitational 


effects of the Sun 



55 


SUBROUTINE DRAG1: Calculates changes In orbital elements due to atmospheric 

drag. 

FUNCTION RHA: Calculates atmospheric density (used in atmospheric drag 

calculations). 

SUBROUTINE CHANGE: Sums all changes in orbital -lements and calculates new 

orbit parameter. 

SUBROUTINE JULCAL: Converts Julian date to calendar date. 

SUBROUTINE CALJUL: Converts calendar date to Julian date. 

SUBROUTIN'. uUXIN: Reads solar flux and geomagnetic index predictions from 

flux file. 


i 
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APPENDIX B. SOLAR FLUX OATA 

Table B.l is a copy of the flux data file created from the June 1985 
update of solar flux and geomagnetic index data obtained from the NASA 
Marshall Space Flight Center (MSFC). The first line consists of seven 
numbers: the year and month the data begin (January 1984), the year and month 

the data end (November 1998), the year and month the data were updated by MSFC 
(June 1985), and the number of data points (179). Each of the remaining lines 
contains the year and month (day is assumed to be the 15th), followed by the 
+2 a and nominal predictions of the 10.7 cm solar flux F^ ^ and the +2 <? and 
nominal predictions of the geomagnetic index, Ap, respectively. In the table, 
for the dates January 1984 through May 1985, actual measured monthly means of 
F io , 7 and Ap appear. Because these are actual values and not predictions, 
there is no distinction between nominal and +2 a values. The predictions begin 


for June 1985. 



TABLE B. 1 . - FLUX DATA FILE 
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APPENDIX C. ORBITAL INPUT OPTIONS 

Symbols 

Az Inertial azimuth, measured positive clockwise from north, deg. 
a Semimajor axis of orbit, R e 
e Eccentricity of orbit, dimensionless 

h fl Height of apogee above Earth, km 
hp Height of perigee above Earth, km 
i Inclination of orbit, deg. 

K OL program FORTRAN array containing integer data 
p Semilatus rectum of orbit, R e 

R fi Earth radius, 6.378145 x 10^ cm 

R.j Initial geocentric radius, km 

v True anomaly, deg. 

Inertial velocity of spacecraft in direction of flight, km/sec 
V Q Square of the fraction of local circular orbit velocity at R^, 
dimensionless 

y Inertial flight path angle, deg. 

aTq Number of days since January 1.0, 1965 Universal Time, days 
At Fr t ^ion of day measured from Greenwich midnight, days 
a Initial longitude, positive east of Greenwich, deg. 
x. Longitudinal increment between ascending node and initial orbit point, 
deg. 

O O 

y Product of universal gravitational constant and Earth mass, cur /sec 
<f> Initial latitude, deg. 

0 Right ascension of ascending node, deg. 

Instantaneous right ascension of Greenwich meridian, deg. 
w Argument of perigee, deg. 
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Calculation of Orbital Elements from Input 

As mentioned in the Input Format section, there are three ways to input 
the orbit parameters, corresponding to three values of K(50). 

K(50) = 0: Input the orbital elements p, e, o) , ft 9 and 1 

K(50) = 1: Input u, i, X, h and h. 

K{50) = -I: Input the trajectory variables Az, <j,, x, R^ V^, and y. 

If K(50) * I or K{50) = -I, the orbital elements p, e, ca, n, and i are 
calculated as follows (see fig. C.I for orbit geometry): 

For K(50) = 1 

e: a = (h p + h a )*10 5 /(2R e ) 

e * (h a *10 5 /R e + l)/a - 1 

p: p » a(l-e 2 ) 

w: Input 

P: = 100.43735 + 0.98561 aT_ + 360.98561 At 

X i = cos ” 1 (cos (to + v)/cos <p ) 
n X x^. 

i : Input 


For K(50) » -1 

e: V Q = (V i *10 5 ) 2 (R i *10 5 )/ y 

e = (1-(2-V 0 )V 0 cos 2 y ) 1/2 

P : a = -y /(((V.* 105)2 . 2 u /(R.n0S))*R e ) 

p = a(l-e 2 ) 

uj: u - sin'^sin 4 >/sin i)-v 

n: Same as for K(50) * 1 

i : i = cos ' 1 (sin Az cos <j») 


— a - 



Gree. wich - 
Meridian 


63 



Figure C.l. - Orbit geometry for input. 
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